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Abstract 


In  recent  work  we  have  formnlated  a  new  approach  to  compressible  flow  simulation, 
combining  the  advantages  of  shock-fitting  and  shock-capturing.  Using  a  cell-centered 
Roe  scheme  discretization  on  unstructured  meshes,  we  warp  the  mesh  while  marching 
to  steady  state,  so  that  mesh  edges  align  with  shocks  and  other  discontinuities.  This 
new  algorithm,  the  Shock-fitting  Lagrangian  Adaptive  Method  (SLAM)  is,  in  effect,  a 
reliable  shock-capturing  algorithm  which  yields  shock-fitted  accuracy  at  convergence. 

Shock-capturing  algorithms  like  this,  which  warp  the  mesh  to  yield  shock-fitted 
accuracy,  are  new  and  relatively  untried.  However,  their  potential  is  clear.  In  the 
context  of  sonic  booms,  accurate  calculation  of  near-field  sonic  boom  signatures  is 
critical  to  the  design  of  the  High  Speed  Civil  Transport  (HSCT).  SLAM  should  allow 
computation  of  accurate  N-wave  pressure  signatures  on  comparatively  coarse  meshes, 
significantly  enhancing  our  ability  to  design  low-boom  configurations  for  high  speed 


*Tbis  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Con¬ 
tract  No.  NASl-19480,  while  the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center, 
Hampton,  VA  23681. 
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1  Introduction 

One  of  the  principal  difficulties  in  computing  compressible  flows  is  that  such  flows  are 
generally  only  piecewise  smooth.  The  solutions  are  smooth,  except  along  a  sequence  of 
arcs  or  surfaces  at  which  the  solution  or  its  derivatives  have  jump  discontinuities.  In  the 
vicinity  of  these  discontinuities,  difference  approximations  are  problematic.  Moreover,  errors 
at  shocks  can  contaminate  the  solution  everywhere. 

There  are  two  basic  approaches  to  computation  of  compressible  flows,  shock-capturing 
and  shock-fitting.  Shock-capturing,  in  which  one  applies  a  well  chosen  difference  scheme 
throughout  the  flow  field,  is  effective  and  reliable,  but  is  usually  only  first  order  accurate 
near  shocks.  Such  schemes  smear  shocks  over  several  mesh  cells,  limiting  the  accuracy  and 
resolution  obtainable. 

The  alternative  is  shock-fitting,  in  which  the  shocks  are  treated  as-  internal  boundaries 
in  the  flow  across  which  one  applies  the  Rankine-Hugoniot  jump  conditions.  Shock-fitting 
algorithms  can  achieve  an  arbitrarily  high  order  of  accuracy,  though  properly  locating  shocks 
is  difficult,  especially  for  flows  containing  complex  embedded  shocks. 

In  recent  work  we  have  formulated  a  new  approach  to  compressible  flow  simulation, 
combining  the  advantages  of  shock-fitting  and  shock-capturing.  The  fundamental  difficulty 
in  shock-fitting  has  always  been  that  of  unambiguously  detecting  and  locating  shocks.  In 
simple  cases,  such  as  that  of  a  strong  bow  shock,  one  has  enough  apriori  knowledge  of 
the  shock  location  that  fitting  schemes  are  highly  successful.  However,  in  more  complex 
situations,  shock-fitting  becomes  difficult  and  unreliable.  For  this  reason,  given  the  simplicity 
and  effectiveness  of  modern  shock- capturing  schemes,  the  latter  have  come  to  dominate 
computational  aerodynamics,  despite  their  limited  resolution. 

The  new  approach  we  are  exploring  begins  with  a  cell-centered  Roe  discretization,  on 
unstructured  meshes  [Ij.  Roe’s  scheme  is  a  popular  and  effective  method,  which  has  an 
interesting  property:  at  steady  state,  this  scheme  imposes  the  exact  Rankine-Hugoniot  jump 
conditions  on  any  cell  face  which  is  oriented  to  and  lying  along  a  shock  or  other  flow  discon¬ 
tinuity.  Thus  if  we  warp  the  mesh  while  marching  to  steady  state,  so  that  shocks  and  other 
discontinuities  lie  along  cell  faces.  Roe’s  scheme  will  give  virtually  exact  answers  there.  This 
is  the  basic  idea  of  the  Shock-fitting  Lagrangian  Adaptive  Method.  SLAM  is,  in  effect,  a 
reliable  shock-capturing  algorithm,  which  incidentally  yields  shock-fitted  solutions  at  conver¬ 
gence,  with  the  attendant  improvement  in  accuracy  and  resolution.  Our  algorithm  is  closely 
related  to  two  other  methods  presented  recently  [2],  [3].  We  give  the  precise  similarities  and 
differences  in  the  next  section. 


2  Related  Work 

While  shock-fitting  has  existed  for  decades  [4,  5,  6],  the  idea  of  warping  an  unstructured 
grid  in  a  shock-capturing  code  to  effectively  fit  shocks  is  new.  The  basic  idea  of  using 
a  conservative  shock-capturing  scheme  on  unstructured  meshes,  and  warping  the  mesh  to 
fit  shocks  during  iteration  to  convergence,  was  independently  developed  by  several  groups, 
including  ourselves. 

All  three  of  the  groups  we  are  aware  of  used  algorithms  based  on  Roe-scheme  discretiza- 
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tions,  but  beyond  that  the  details  of  these  approaches  differ.  Parpia  and  Parikh  [3]  use  the 
waves  occuring  in  a  six-wave  multidimensional  Riemann  solver  to  align  mesh  edges  with 
shocks,  without  actually  fitting  shocks.  The  multidimensional  Riemann  solver,  due  origi¬ 
nally  to  Roe,  is  described  in  references  [7,  8].  Aligning  the  grid  allows  them  to  achieve  true 
“one-point”  shocks,  free  of  the  “splitting  error”  that  occurs  when  shocks  cross  the  mesh  at 
an  oblique  angle.  The  other  group,  Trepanier  et  al.  also  used  the  six-wave  multidimensional 
Riemann  solver  to  control  mesh  warping  [2,  9].  However,  unlike  Parpia  and  Parikh,  they 
also  move  mesh  edges  to  coincide  with  shocks,  thus  obtaining  shock-fitting  accuracy  in  the 
final  solution.  ~ 

Our  algorithm  is  similar  to  that  of  Trepanier  et  al.,  differing  in  that  we  warp  the  mesh 
using  only  density  gradients,  rather  than  using  the  waves  occuring  in  a  multidimensional 
Riemann  solver.  Thus  unlike  these  other  groups,  we  do  not  need  a  separate  discretization 
to  control  mesh  movement.  In  effect,  we  are  reusing  information  from  the  Roe  scheme 
discretization.  Our  approach  is  clearly  simpler  and  cheaper.  Also,  extension  to  3D  is  straight 
forward  with  our  approach.  The  other  approaches  should  extend  to  3D  as  well,  though  how 
well  they  extend  is  not  clear  multidimensional  Riemann  solvers  are  not  nearly  as  robust  and 
effective  in  3D  as  in  2D.  Design  of  better  3D  multidimensional  Riemann  solvers  is  an  area 
of  active  research  [10]. 


3  Algorithm  Design 

The  discretization  used  here  is  the  cell  centered  Roe  scheme.  This  is  an  effective  and 
heavily  used  scheme,  whose  occasional  failings  are  now  well  understood  and  easily  overcome 
[11].  In  our  algorithm  we  march  to  steady  state  using  the  Roe  scheme  coupled  to  a  “locally 
implicit”  time  stepping  scheme  [1].  For  the  first  30  or  40  iterations,  we  keep  the  mesh  frozen, 
allowing  initial  transients  to  dissipate.  After  that,  we  allow  the  mesh  to  warp  at  each  time 
step  to  fit  the  developing  shocks. 


Scheme  A 


Several  schemes  for  warping  the  mesh  have  been  tried.  The  two  found  to  work  best  are 
described  here.  The  first  of  these,  which  we  will  call  “scheme  A”  for  convenience,  is  based 
on  reconstruction  of  ID  shocks  on  the  edges  of  the  mesh.  We  begin  by  interpolating  the 
cell-centered  densities  to  edges  and  vertices.  Edge  values  are  obtained  by  simple  averaging 
of  adjacent  cell  values,  while  vertex  values  are  obtained  by  weighting  the  density  value  in 
each  cell  by  the  angle  subtended  by  that  cell.  That  is,  if  triangles  touch  vertex  v,  the 

angle-weighted  density  at  v  is 


Pv  — 


Et  ^t,v  Pt 
Et  h,. 


where  dT,v  is  the  angle  subtended  by  triangle  T  at  vertex  v.  The  denominator  here  is  360 
degrees,  except  at  boundary  points. 

“Angle-weighting”  is  crucial  to  our  scheme,  since  it  means  that  vertex  values  along  a 
fitted  straight  shock  separating  states  Ui  and  Ur  will  be  correct,  independent  of  the  local 
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mesh  topology  -  Ui  and  Ur  will  both  be  weighted  by  180  degrees.  While  angle-weighting  does 
produce  second  order  errors  along  curved  shocks,  this  has  not  been  found  to  be  a  problem. 

Now  given  density  values  on  all  edges  and  vertices,  we  reconstruct  a  shock  on  every  edge. 
Assume  the  density  profile  along  each  edge  is  piecewise  constant,  with  a  single  jump.  At 
the  end  of  each  edge,  the  density  must  match  that  of  the  vertex,  while  the  jump  location 
is  determined  by  requiring  that  the  average  of  the  reconstructed  piecewise  constant  density 
match  the  known  value  on  this  edge,  if  possible,  as  shown  in  Figure  1. 


Figure  1:  Reconstruction  of  a  Shock  Along  an  Edge 

Once  we  have  reconstructed  shocks  on  the  edges  of  the  mesh,  we  attempt  to  attract  each 
vertex  to  the  shocks  on  the  edges  incident  on  that  vertex.  To  do  this,  we  compute  the  center 
of  mass  of  these  shocks,  assuming  the  mass  of  each  shock  is  proportional  to  its  strength. 
This  center  of  mass  is  a  point  in  the  plane,  to  which  we  attempt  to  attract  this  vertex.  Note 
that  we  are  attracting  only  vertices  to  shocks.  There  is  no  direct  consideration  of  either 
attracting  or  orienting  edges.  In  effect,  we  are  making  use  of  the  following  principle: 

When  two  vertices  of  a  triangle  (or  three  of  a  tetrahedron)  lie  on  a  straight 
shock)  the  intervening  cell  face  exactly  fits  the  shock. 

Scheme  B 

Scheme  A  is  effective,  but  does  not  work  well  for  very  weak  shocks.  Scheme  B  remedies 
this.  It  is  composed  of  two  separate  components: 

1.  A  shock  detector 

2.  A  vertex  attraction  force 

The  latter  is  used  only  at  points  detected  as  shocks. 

The  shock  detector  uses  density  gradients  at  the  vertices,  computed,  for  example,  by 
Green’s  theorem  path  integrals.  Let  g^  denote  the  density  gradient  at  vertex  v  and  time 
step  n.  For  all  neighboring  vertices,  a  of  v,  define  a  weight 

<  =  (a  -  u)  •  < 

where  •  denotes  the  normal  inner  product.  Then  we  take  the  weighted  average  of  gradient 
norms,  with  respect  to  these  weights: 

n  _  'fhneighhoTS  '^a  1  II 

^neighbors 
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Figure  2:  Attraction  of  Vertices  to  Shock 

We  flag  points  at  which  the  density  gradient  exceeds  this  average  of  gradients  at  sur¬ 
rounding  points.  In  particular,  we  threshold  the  quantity 

max{0,  llg^ll  -  c^) 

Ikll  +  ^  ’ 

where  e  >  0  is  needed  to  avoid  division  by  zero  in  smooth  regions.  With  both  numerator 
and  denominator  proportional  to  the  density  gradient  (neglecting  epsilon),  this  detector  is 
equally  effective  at  detecting  weak  and  strong  shocks. 

The  vertex  movement  force  at  vertex  v  in  this  scheme  is  of  the  form: 

^  lk"l! 

That  is,  a  scalar  multiple,  of  the  unit  vector  in  the  direction  of  the  density  gradient.  The 
scale  factor  s": 

_  in  /  -nx  cn 

\Pv  Pv) 

with: 

/i”  .  local  mesh  size 

p"  .  angle- weighted  vertex  average 

.  local  ambient  density 

.  magnitude  of  local  density  range 

We  take  /i"  to  be  the  minimum  length  of  edges  incident  on  v,  while  p"  is  the  angle- weighted 
vertex  average,  as  before. 

Quantities  and  both  depend  on  the  range  of  density  in  a  small  surrounding  region. 
Define  the  “local  maximum  density”  p”^,  as  the  maximum  density  in  cells  touching  vertices 
neighboring  v.  Thus  p"  „  is  the  maximum  density  in  the  20  or  so  surrounding  cells.  Similarly, 
define  the  “local  minimum  density”  pg^^  as  the  minimum  density  in  these  surrounding  cells. 
Then  the  local  density  range  is: 

cn  ^n  jn 

~  Pl,v  Po,v 

Similarly,  the  local  “ambient  density”  is 

p:  =  (pi.  +  piJ/i 
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Note  that  this  is  an  average  of  density  extremes,  as  opposed  to  a  direct  average  of  densities. 

The  motivation  behind  all  of  this  is  the  following.  A  vertex  along  a  shock  is  correctly 
located  when  the  density  value  there  is  midway  the  high  and  low  density  in  a  surrounding 
region.  Thus  we  want  to  satisfy  the  equation 


at  vertices  along  shocks.  The  scale  factor  s”  approximates  the  amount  of  the  correction 
needed  to  satisfy  this  equation.  Since  we  adjust  the  grid  at  every  iteration,  the  precise  scale 
factor  used  is  unimportant;  in  effect  it  is  just  a  relaxation  parameter. 


Mesh  Control 

Using  either  scheme  A  or  scheme,  vertices  are  rapidly  attracted  to  shocks.  However, 
without  constraints  on  mesh  movement,  one  rapidly  produces  undesirably  thin  cells,  or 
negative  cell  areas.  To  avoid  this  two  things  are  needed: 

1.  mesh  control  forces,  partially  counteracting  the  forces  attracting  shocks  to  vertices. 

2.  control  of  mesh  movement  step-size  control. 

We  are  currently  using  two  forces,  one  proportional  to  the  change  in  cell  area,  another 
based  on  angles  at  vertices.  The  latter,  which  applies  torques  on  edges,  prevents  angles 
from  approaching  either  0  or  180  degrees.  As  angles  approach  either  0  or  180  degrees,  the 
torques  it  produces  become  infinite.  Thus  if  the  ODE  governing  mesh  movement  is  properly 
integrated,  degenerate  triangles  cannot  occur.  By  contrast,  using  “springs”  on  edges  is  not 
effective,  since  they  cannot  prevent  angles  from  approaching  0  or  180  degrees. 

There  are  a  number  of  ways  to  control  the  step-size  in  the  mesh  movement  scheme.  Our 
approach  is  to  compute  a  maximum  step  at  each  vertex,  designed  to  prevent  degenerate 
triangles.  For  every  triangular  cell  R  let  hminiT)  be  the  minimum  of  the  three  altitudes. 
Then  for  each  vertex  define 

~  min 

neighboring  triangles 

If  no  vertex  v  moves  further  than  ^hmin{v)i  degenerate  triangles  cannot  occur. 

Note  that  controlling  step-size  alone  suffices  to  prevent  degenerate  triangles,  but  grid 
quality  may  be  quite  poor.  The  combination  of  these  mesh  control  forces  and  step-size 
control  suffices  to  maintaining  mesh  quality,  while  still  allowing  effective  fitting  of  shocks. 

Generalized  Van  Albada  Limiter 

The  first  order  scheme  just  described  works  well,  but  provides  inadequate  resolution  in 
smooth  regions.  Second  or  third  order  accuracy  can  be  achieved  with  a  MUSCL-style  scheme 
[12],  in  which  one  reconstructs  a  polynomial  in  every  cell  via  an  appropriate  “limiter.”  One 
way  of  doing  this  is  to  adapt  the  stencils,  following  the  ENO  approach.  However  ENO  is 
complex  and  expensive  on  unstructured  meshes  [13]. 
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Our  approach  is,  instead,  to  use  a  multidimensional  generalization  of  the  Van  Albada 
limiter  [14].  This  limiter  is  simple,  reliable,  and  has  the  attractive  property  of  not  clipping 
extrema.  Thus  it  can,  in  principle,  achieve  perfectly  sharp  approximations  of  N-waves  on 
very  coarse  meshes. 

The  goal  in  a  MUSCL  scheme  is  to  replace  the  constant  value  in  each  cell  by  a  linearly 
varying  distribution 

=  qf  +  i‘n-Vi){^q)7 

where  is'an  approximation  to  the  gradient.  The  Van  Albada  limiter  takes  this  gradient 
as  a  nonlinear  average  of  the  gradients  computed  by  forward  and  backward  differencing: 


using  the  averaging  function 


ave  (a,  b) 


(6^  +  e^)  a  +  (a^  +  e^)  b 
+  62  2  £2  ’ 


where  epsilon  is  a  small  positive  constant  designed  to  provide  smooth  transitions.  This  kind 
of  averaging  was  used  in  [14]  for  all  quantities  except  density,  which  was  handled  slightly 
differently  to  avoid  negative  overshoots  in  strong  astrophysical  flows. 

The  Van  Albada  limiter  generalizes  easily  to  unstructured  meshes.  To  see  this,  rewrite 
the  above  formulas  as 


{Sq)^  =  ^a(g",i-9r)  +  y^b{q?-qlr), 

with 

{P  + 

“  a2  +  62  +  2  £2’ 

(a^  +  P) 

a2  +  62  +  2  £2 

Now  in  a  similar  way,  for  triangular  mesh  cells,  assume  one  has  gradients  g",  g" 
vertices  of  a  triangle,  obtained,  as  usual,  by  Green’s  theorem  path  integrals.  Then  one  can 
compute  the  cell  centered  gradient  as 

g"  =  g"  +  Wb  gb  +  Wo  g" 

for  suitable  weights  u>a,  rut,  Wc-  Constraining  these  weights  by 

Wa  +  Wb  +  Wc  =  1 

Wa,  Wb,  Wc  e  [0,1]. 

yields  second  order  consistency  of  the  overall  scheme,  cissuming  the  nodal  gradients  are  first 
order  accurate.  We  also  want  to  preserve  the  Van  Albada  property  of  not  clipping  extrema. 
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The  particular  choice  we  used  was 


(h 

c 

+ 

a  b 

+ 

b 

c 

+ 

c  a 

+  3  £2 

(c 

a 

+ 

a  b 

+ 

b 

c 

+ 

c  a 

-f3  £2 

{a 

b 

+ 

a  b  +  be  +  ca  +  3e^’ 


with  a  =  i  c  —  IkelP-  Other  choices  work  about  as  well.  In  particular, 

one  can  chose 

«  =  Ikll,  b  =  il^fcll,  c  =  ll^cll, 

in  closer  analogy  with  the  original  Van  Albada  scheme.  We  prefer  the  stronger  switching 
that  occurs  in  using  the  squares. 

This  generalized  Van  Albada  limiter  has  the  property  that  near  strong  jumps  the  re¬ 
constructed  gradients  use  information  entirely  from  one  side  of  the  jump,  thus  achieving 
second-order  consistency  while  avoiding  spurious  oscillations.  Thus  one  can  think  of  this  as 
an  inexpensive  approach  to  ENO,  avoiding  the  use  of  complex  adaptive  stencils  and  to  some 
extent  the  convergence  difficulties  they  create. 

Parabolic  Reconstruction 

An  interesting  question  to  ask  is  whether  the  same  approach  can  be  used  for  higher  order 
reconstruction.  While  we  have  implemented  only  the  second,  it  is  clear  one  can  easily  go  to 
third  or  higher  order.  Given  the  original  piecewise  constant  states  in  each  cell,  the  Green’s 
theorem  path  integrals  yield  first  order  accurate  gradients  at  the  vertices,  so  the  resulting 
piecewise  linear  state  in  each  cell  is  second  order  accurate.  We  can  at  this  point  repeat 
the  path  integrals  to  evaluate  the  gradients  at  the  vertices  to  second  order,  or  the  second 
derivatives  at  the  vertices  to  first  order.  The  resulting  information  more  than  suffices  for 
a  Van  Albada-style  parabolic  reconstruction  in  each  cell,  of  third  order  accuracy,  and  the 
process  can  be  continued. 

Though  this  “boot  strapping”  approach  might  at  first  seem  suspect,  it  is  perfectly  valid. 
In  repeating  the  path  integrations,  one  is  effectively  widening  the  difference  stencil,  thus 
allowing  higher  order  approximation.  There  are,  however,  a  couple  of  difficulties  here.  First, 
the  higher  order  path  integrals  required  are  more  expensive.  Second,  to  fully  exploit  higher 
order  reconstructions,  one  has  to  compute  the  flux  balances  via  higher  order  Gauss  formu¬ 
las.  For  example,  with  parabolic  reconstruction,  one  needs  two  quadrature  points  per  side, 
doubling  the  cost  of  the  method.  However,  even  without  higher  order  quadrature  of  the 
fluxes,  it  may  still  make  sense  to  perform  parabolic  reconstruction,  as  in  the  “Piecewise 
Parabolic  Method”  of  Woodward  and  Collela  [15]  which,  in  effect,  reduces  the  constant  in 
the  truncation  error  while  remaining  formally  second  order  accurate. 
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4  Sonic  Boom  Computations 

While  SLAM  can  be  applied  to  steady-state  Euler  calculation  (and  presumably  steady- 
state  Navier  Stokes  calculations  as  well)  one  application  for  which  it  is  especially  well  suited 
is  computation  of  sonic  booms.  Any  aircraft  moving  faster  than  the  speed  of  sound  generates 
a  sonic  boom.  Nefir  the  aircraft,  the  various  shocks,  rarefactions,  wakes  and  boundary  layers 
interact  in  complex,  nonlinear  ways.  However,  a  few  body-lengths  from  the  aircraft,  these 
nonlinear  waves  usually  coalesce  into  a  characteristic  N-wave,  which  then  propagates  linearly 
to  the  ground.  As  it  does  so,  this  N-wave  undergoes  diffraction  through  the  atmospheres 
density  gradients,  and  Doppler  shifts  in  the  wind  layers  encountered. 

One  can  compute  the  shock  waves  near  the  aircraft  with  any  three  dimensional  flow  solver. 
However,  current  mesh  resolutions  are  inadequate  for  detailed  sonic  boom  calculation  [16, 17]. 
Moreover,  shock  steepening  is  still  occurring  many  body-lengths  from -the  aircraft,  so  that 
a  very  large  mesh  may  be  necessary.  Through  Lagrangian  mesh  movement,  SLAM  can  find 
even  extremely  weak  shocks  and  compute  them  accurately  on  relatively  coarse  meshs.  Thus 
its  potential  capabilities  in  this  kind  of  computation  should  be  much  better  than  competing 
adaptive-mesh  methods. 

Though  SLAM  is  not  yet  implemented  in  3D,  its  theoretical  potential  is  easy  to  assess. 
One  can  make  the  argument  for  the  dramatic  increase  in  efficiency  of  Lagrangian  fitting 
strategies  over  mesh-enrichment  strategies  through  a  simple  geometric  argument,  like  that 
given  by  Swartz  [18].  Instead,  we  look  here  at  a  simple  model  problem.  Take  as  model 
problem  the  problem  of  calculating  the  bow  shock  from  a  high  speed  projectile.  Assume,  for 
convenience  the  shock  is  a  perfect  cone,  and  assume  we  wish  to  resolve  the  shock  in  cells  of 
size  6  in.  a  domain  of  radius  R,  as  shown; 


Figure  3:  Resolving  a  Conical  Shock 


Using  a  uniform  mesh,  the  number  of  cells  required,  N,  evidently  satisfies 


N  oc 


or  equivalently: 


S  oc 
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Adaptive  refinement  strategies  greatly  improve  the  resolution  achieved  with  a  given  num¬ 
ber  of  cells.  Consider,  for  example,  a  Cartesian  adaptive  strategy  based  on  oct-tree  refine¬ 
ment.  In  this  case,  the  total  number  of  cells  required  is  proportional  to  the  number  required 
to  tile  the  shock  (assuming  cells  sizes  grow  geometrically  away  from  the  shock).  Thus  we 
have 

-  -  (f) 


S  oc  - 

In  effect  the  problem  dimension  has  been  reduced  by  one. 

With  SLAM,  one  can  resolve  the  shock  to  within  a  distance  of  S  using  tetrahedrons  of 
much  larger  diameter.  The  positional  error  in  approximating  a  fitted  shock  with  tetrahedrons 
of  diameter  h  is 

e  K 

where  k  is  the  local  curvature  of  the  shock.  Then  setting  e  =  6  we  find  we  should  take: 


h  =  \f^T^ 


Now  the  curvature  /c  of  a  conical  shock  goes  as  1/r,  where  r  is  the  distance  from  the  axis. 
Thus  the  size  of  cells  in  a  properly  graded  mesh  would  be: 

h{r)  = 

Now  to  count  the  cells  tiling  the  shock,  define  the  cell  density  d{r)  =  lfh^{r).  The 
number  of  tetrahedrons  needed  to  tile  the  shock  is  then  the  integral  of  d{r)  over  the  shock 

fR  rR 

N  (X  2-Krdir)  dr  =  27r  -dr 
Jo  Jo  0 


iV  oc 


6  oc  RN- 


The  total  number  of  cells  will  again  be  proportional  to  the  number  required  to  tile  the  shock. 
Thus  we  have  the  following  surprising  result: 


The  cost  of  resolving  a  conical  shock  wave  in  three  dimensions  to  within  6  using 
SLAM  is  asymptotically  the  same  as  that  of  solving  a  one  dimensional  shock 
problem  on  a  uniform  mesh  of  spacing  S. 


In  other  words,  the  combination  of  fitting  together  with  mesh-grading  is  sufficient  to  achieve 
third  order  global  accuracy,  despite  the  point-singularity  at  the  nose. 

This  is  a  theoretical,  asymptotic  result,  neglecting  the  many  complex  implementation 
issues  involved.  It  also  assumes  that  the  flow  near  the  body  can  be  solved  with  cost  at  most 
proportional  to  the  cost  of  computing  the  shock,  and  that  the  smooth  regions  in  the  far  field 
can  be  resolved  on  such  coarse  meshes  that  their  cost  is  negligible.  Modulo  these  caveats, 
this  analysis  is  valid,  and  suggests  that  computation  of  sonic  booms  produced  by  supersonic 
aircraft  should  be  quite  tractable,  with  SLAM-style  algorithms. 
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5  Experimental  Results 

Results  obtained  with  SLAM,  while  preliminary,  are  very  promising.  We  have,  in  general, 
no  trouble  with  strong  shocks,  including  attached  and  detached  bow  shocks,  fish  tail  shocks, 
and  standing  shocks  on  transonic  airfoils.  Similar  techniques  can  be  used  to  resolve  slip  lines 
and  contacts  [9],  though  we  have  not  yet  studied  this. 

Figure  4  shows  an  un-adapted  mesh  of  8,000  points  around  a  10%  circular  arc  airfoil, 
while  Figure  5_ shows  the  same  mesh  after  modification  by  SLAM  to  align  mesh  edges  with 
shocks.  Figure  6  shows  the  density  field  on  this  8,000  point  adapted  mesh.  The  shocks  are 
sharp  all  the  way  to  the  far-field  boundary,  and  the  limiter  is  also  producing  an  accurate 
solution  in  the  smooth  region  between  shocks. 

One  can  judge  the  solution  more  accurately  by  taking  cross  sections.  Figures  [7-10]  show 
SLAM  solutions  on  the  8,000  point  mesh  in  Figure  5,  and  on  an  analogous  4,000  point  mesh. 
Figures  7  and  8  show  the  density  cross  sections  taken  parallel  to  the  axis  one  chord  from 
the  axis,  while  Figures  9  and  10  show  density  cross  sections  four  chords  from  the  axis.  As 
can  be  seen,  the  sonic-boom  profile  is  well  captured,  even  four  chords  from  the  axis.  These 
solutions  are  mesh-converged  to  graphical  accuracy,  in  smooth  regions.  However,  there  is  a 
slight  anticipation  of  the  bow  shock,  worse  on  the  4,000  point  mesh,  due  to  error  in  shock 
location.  There  are  several  numerical  effects,  creating  second-order  errors  in  shock  location. 

For  this  simple  problem,  one  can  obtain  qualitatively  reasonable  solutions  via  SLAM, 
using  fewer  than  1,000  mesh  points.  However,  accuracy  is  lacking  until  the  smooth  regions 
are  resolved.  By  contrast.  Figures  [11-14]  show  the  solution  on  the  un-adapted  8,000  point 
mesh  of  Figure  4,  and  on  an  analogous  32,000  points  mesh.  With  the  Lagrangian  mesh 
adaptivity  turned  off,  the  sonic  boom  profiles  are  now  badly  distorted,  especially  on  the 
8,000  point  mesh,  even  though  we  made  a  real  effort  to  concentrate  mesh  in  regions  where 
the  sonic  boom  was  expected. 

On  the  8,000  point  mesh,  at  four  chords  from  the  axis.  Figure  13,  the  bow  and  tail 
shocks  are  separated  by  about  15  mesh  widths.  Thus  significant  smearing  is  inevitable.  This 
smearing  is  not  as  severe  on  the  32,000  point  mesh,  which  has  four  times  the  mesh  density 
throughout  the  flow  field,  but  the  answer  there  is  still  much  poorer  than  the  SLAM  solutions. 
In  particular,  comparing  Figures  14  and  10,  notice  that  the  extrema  are  substantially 
blurred  on  the  32,000  point  un-adapted  mesh.  The  combination  of  Lagrangian  adaptivity 
and  generalized  Van  Albada  limiter  is  particularly  effective  at  getting  extrema  right. 

For  this  simple  example,  one  requires  at  least  60  times  more  points  with  an  un-adapted 
mesh  than  with  SLAM  for  comparable.  For  example,  halving  the  mesh  size  in  the  32,000 
point  mesh,  one  would  yield  a  128,000  point  mesh  achieving  accuracy  approaching  that  of 
the  4,000  point  SLAM  solution. 

Figure  15  shows  a  somewhat  more  complicated  example,  flow  over  an  airfoil  with  a 
perfectly  sharp  nose.  Since  the  interior  angle  at  the  nose  of  this  airfoil  is  zero,  there  is  no 
shock  there.  Instead,  a  shock  forms  in  the  free  stream,  some  distance  away,  where  the  acoustic 
waves  coalesce.  Figures  16  and  17  show  the  cross  section  0.1  and  0.2  chords  respectively 
from  the  axis.  The  smooth  profile  at  the  nose  in  Figure  16  has  begun  to  sharpen  into  a 
shock  in  Figure  17.  Steepening  continues,  as  shown  in  Figure  18,  at  0.4  chords  from  the 
axis.  At  0.8  chords,  Figure  19,  the  eventual  N-wave,  is  well  established. 

The  point  here  is  that,  unlike  most  shock-fitting  schemes,  the  Lagrangian  adaptive  scheme 
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has  no  effect  on  the  underlying  conservative  discretization.  Thus  examples  like  this,  with 
coalescing  waves,  intersecting  shocks,  and  so  on,  present  no  difficulty,  at  least  in  principal. 

We  have  not  yet  done  experiments  to  compare  SLAM  with  a  standard  mesh-enrichment 
adaptive  strategy,  though  there  is  little  doubt  how  the  comparison  will  come  out.  SLAM 
achieves  shock-fitted  accuracy  without  addition  of  mesh  points,  while  enrichment  strategies 
substantially  increase  in  the  number  of  mesh  points.  Moreover,  since  enrichement  strategies 
leave  the  mesh  unaligned,  the  well  resolved  shocks  there  will  still  be  smeared  over  several 
mesh  points.  Thus  SLAM  should  require  about  one  tenth  as  many  mesh  points  as  mesh- 
enrichment  schemes  in  2D,  and  should  be  relatively  even  better  in  3D.  The  results  of  such  a 
comparison  will  be  reported  in  a  sequel. 

6  Conclusions 

Lagrangian  adaptive  grid  methods,  like  SLAM,  have  great  potential  for  resolving  shocks 
and  other  flow  discontinuities.  Unlike  mesh-enrichment  strategies,  which  put  much  finer  mesh 
along  shocks,  fitting  strategies  can  resolve  discontinuities  without  increasing  the  number  of 
mesh  cells.  This  advantage  is  especially  important  in  three  dimensions,  where  the  cost  of 
tiling  shocks  with  fine  mesh  is  great.  This  improved  resolution  is  achieved  at  little  cost, 
and  without  loss  of  robustness,  since  we  retain  a  Roe  scheme-based  shock  capturing  scheme. 
Thus  even  if  the  fitting  scheme  fails,  we  still  have  a  robust  and  effective  shock-capturing 
algorithm. 

We  have  demonstrated  the  efficacy  of  SLAM  in  2D,  and  are  beginning  work  on  an  anal¬ 
ogous  3D  code.  The  latter  is  intended  to  be  applied  to  the  problem  of  predicting  the  sonic 
boom  profiles  of  supersonic  aircraft.  Current  CFD  codes  cannot  adequately  resolve  the  com¬ 
plex  shock  waves  emanating  from  a  supersonic  vehicle,  since  one  cannot  afford  a  sufficiently 
fine  grid  extending  4-5  body-lengths  from  the  aircraft.  Fitting  schemes,  like  SLAM,  will 
be  able  to  do  much  better,  and  should  be  able  to  reproduce  the  complex  shock  patterns 
observed  in  wind-tunnel  experiments. 
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Figure  6:  Flow  Field,  8000  Point  Adapted  Mesh 
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Figure  7:  Density  One  Chord  from  Axis,  Circular  Arc  Airfoil 
Mach  1.4,  4000-point  Adapted  Mesh 


Figure  8:  Density  One  Chord  from  Axis,  Circular  Arc  Airfoil, 
Mach  1.4,  8000-point  Adapted  Mesh 
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Figure  9:  Density  Four  Chords  from  Axis,  Circular  Arc  Airfoil 
Mach  1.4,  4000-point  Adapted  Mesh 
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Figure  10:  Density  Four  Chords  from  Axis,  Circular  Arc  Airfoil 
Mach  1.4,  8000-point  Adapted  Mesh 
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Figure  13:  Density  Four  Chords  from  Axis,  Circular  Arc  Airfoil 
Mach  1.4,  8000-point  Un-adapted  Mesh 
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Figure  14:  Density  Four  Chords  from  Axis,  Circular  Arc  Airfoil 
Mach  1.4,  32000-point  Un-adapted  Mesh 
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Figure  16:  Density  0.1  Chords  from  Axis,  Pointed-nose  Airfoil 
Mach  1.8,  9000-point  Adapted  Mesh 


Figure  17:  Density  0.2  Chords  from  Axis,  Pointed-nose  Airfoil 
Mach  1.8,  9000-point  Adapted  Mesh 
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Figure  18:  Density  0.4  Chords  from  Axis,  Pointed-nose  Airfoil 
Mach  1.8,  9000-point  Adapted  Mesh 
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Figure  19:  Density  0.8  Chords  from  Axis,  Pointed-nose  Airfoil 
Mach  1.8,  9000-point  Adapted  Mesh 
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